##### R CODE TO REPLICATE THE FIGURES AND TABLES INCLUDED IN 
##### Adrian Lucardi, "The Effect of District Magnitude on Electoral Outcomes. Evidence from Two Natural Experiments in Argentina," British Journal of Political Science, 49(2), 2019, 557-577.

## Loading the required R packages. Make sure to have them installed before proceeding. In order to install a package called X, write install.packages("X")
library (car)
library (clusterSEs)
library (foreign)
library (ggplot2)
library (gridExtra)
library (lattice)
library (lmtest)
library (plm)
library (plyr)
library (sandwich)
library (scales)
library (survival)
library (xtable)

## display options
options (digits=4, scipen=8, show.signif.stars=FALSE)

## working directory -> change this to the directory where your have your data
setwd ("/Users/adrianlucardi/Dropbox/Replication files/")



###### DOWNLOADING AND PROCESSING THE DATA

## Downloading the data
baseFull <- read.csv ("Argentine legislative elections 1983-2015.csv", sep=",", header=TRUE)
baseFull$year2 <- with (baseFull, factor (paste (year, ffCode, sep="")))  ## we need this b/c the counterfactuals have the same year
summary (baseFull); dim (baseFull)


## Creating the .plm datasets --> needed for the FE models

# actual values
base <- pdata.frame (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], index=c ("district", "year"))

# small provinces (delegation < 5)
base.sm <- pdata.frame (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], index=c ("district", "year"))

## counterfactuals to calculate the mechanical effect (A -> B)
baseAB <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="A" | baseFull$ffCode=="B") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], index=c ("district", "year2"))

# small provinces (delegation < 5)
baseAB.sm <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="A" | baseFull$ffCode=="B") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], index=c ("district", "year2"))

# counterfactuals to calculate the psychological effect (B -> D)
baseBD <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="B" | baseFull$ffCode=="D") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], index=c ("district", "year2"))

# small provinces (delegation < 5)
baseBD.sm <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="B" | baseFull$ffCode=="D") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], index=c ("district", "year2"))

# counterfactuals to calculate the components psychological effect (I) (A -> C)
baseAC <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="A" | baseFull$ffCode=="C") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], index=c ("district", "year2"))

# small provinces (delegation < 5)
baseAC.sm <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="A" | baseFull$ffCode=="C") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], index=c ("district", "year2"))

# counterfactuals to calculate the components psychological effect (II) (C -> D)
baseCD <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="C" | baseFull$ffCode=="D") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], index=c ("district", "year2"))

# small provinces (delegation < 5)
baseCD.sm <- pdata.frame (baseFull[!is.na (baseFull$ffCode) & (baseFull$ffCode=="C" | baseFull$ffCode=="D") & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], index=c ("district", "year2"))



########## FITTING THE MODELS

## setting the seed and the number of replications for the bootstraps
set.seed (384292)
B <- 999


### Replicating the models reported in Table A2 and Figure 3

## Panel (a): Main effect

# DV is # lists running
(sum.m211a.fe <- summary (m211a.fe <- plm (nLists ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m211a.fe <- coeftest (m211a.fe, vcov=vcovHC (m211a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m211a.fe2 <- summary (m211a.fe2 <- plm (nLists ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))  ## same model, using year factor for computational reasons
(ci.m211a.fe2 <- cluster.wild.plm (m211a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m211d.fe <- summary (m211d.fe <- plm (nLists ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m211d.fe <- coeftest (m211d.fe, vcov=vcovHC (m211d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m211d.fe2 <- summary (m211d.fe2 <- plm (nLists ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m211d.fe2 <- cluster.wild.plm (m211d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

# DV is ENPV
(sum.m221a.fe <- summary (m221a.fe <- plm (enpv ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m221a.fe <- coeftest (m221a.fe, vcov=vcovHC (m221a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m221a.fe2 <- summary (m221a.fe2 <- plm (enpv ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))
(ci.m221a.fe2 <- cluster.wild.plm (m221a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m221d.fe <- summary (m221d.fe <- plm (enpv ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m221d.fe <- coeftest (m221d.fe, vcov=vcovHC (m221d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m221d.fe2 <- summary (m221d.fe2 <- plm (enpv ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m221d.fe2 <- cluster.wild.plm (m221d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

# DV is % vote for two largest parties
(sum.m231a.fe <- summary (m231a.fe <- plm (vote02sum ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m231a.fe <- coeftest (m231a.fe, vcov=vcovHC (m231a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m231a.fe2 <- summary (m231a.fe2 <- plm (vote02sum ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))
(ci.m231a.fe2 <- cluster.wild.plm (m231a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m231d.fe <- summary (m231d.fe <- plm (vote02sum ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m231d.fe <- coeftest (m231d.fe, vcov=vcovHC (m231d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m231d.fe2 <- summary (m231d.fe2 <- plm (vote02sum ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m231d.fe2 <- cluster.wild.plm (m231d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))


## Panel (b): Heterogeneous effects

# DV is # lists running
(sum.m212a.fe <- summary (m212a.fe <- plm (nLists ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m212a.fe <- coeftest (m212a.fe, vcov=vcovHC (m212a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m212a.fe <- vcovHC (m212a.fe, method="white2", type="HC3", cluster="group"))

(sum.m212d.fe <- summary (m212d.fe <- plm (nLists ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m212d.fe <- coeftest (m212d.fe, vcov=vcovHC (m212d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m212d.fe <- vcovHC (m212d.fe, method="white2", type="HC3", cluster="group"))

# DV is ENPV
(sum.m222a.fe <- summary (m222a.fe <- plm (enpv ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m222a.fe <- coeftest (m222a.fe, vcov=vcovHC (m222a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m222a.fe <- vcovHC (m222a.fe, method="white2", type="HC3", cluster="group"))

(sum.m222d.fe <- summary (m222d.fe <- plm (enpv ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m222d.fe <- coeftest (m222d.fe, vcov=vcovHC (m222d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m222d.fe <- vcovHC (m222d.fe, method="white2", type="HC3", cluster="group"))

# DV is % vote for two largest parties
(sum.m232a.fe <- summary (m232a.fe <- plm (vote02sum ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m232a.fe <- coeftest (m232a.fe, vcov=vcovHC (m232a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m232a.fe <- vcovHC (m232a.fe, method="white2", type="HC3", cluster="group"))

(sum.m232d.fe <- summary (m232d.fe <- plm (vote02sum ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m232d.fe <- coeftest (m232d.fe, vcov=vcovHC (m232d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m232d.fe <- vcovHC (m232d.fe, method="white2", type="HC3", cluster="group"))


## Panel (c): Rank-based DV

# DV is # lists running
(sum.m213a.fe <- summary (m213a.fe <- plm (nLists.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m213a.fe <- coeftest (m213a.fe, vcov=vcovHC (m213a.fe, method="white2", type="HC3", cluster="group")))

(sum.m213d.fe <- summary (m213d.fe <- plm (nLists.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m213d.fe <- coeftest (m213d.fe, vcov=vcovHC (m213d.fe, method="white2", type="HC3", cluster="group")))

# DV is ENPV
(sum.m223a.fe <- summary (m223a.fe <- plm (enpv.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m223a.fe <- coeftest (m223a.fe, vcov=vcovHC (m223a.fe, method="white2", type="HC3", cluster="group")))

(sum.m223d.fe <- summary (m223d.fe <- plm (enpv.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m223d.fe <- coeftest (m223d.fe, vcov=vcovHC (m223d.fe, method="white2", type="HC3", cluster="group")))

# DV is % vote for two largest parties
(sum.m233a.fe <- summary (m233a.fe <- plm (vote02sum.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m233a.fe <- coeftest (m233a.fe, vcov=vcovHC (m233a.fe, method="white2", type="HC3", cluster="group")))

(sum.m233d.fe <- summary (m233d.fe <- plm (vote02sum.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m233d.fe <- coeftest (m233d.fe, vcov=vcovHC (m233d.fe, method="white2", type="HC3", cluster="group")))


## Panel (d): Heterogeneous effects (rank)

# DV is # lists running
(sum.m214a.fe <- summary (m214a.fe <- plm (nLists.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m214a.fe <- coeftest (m214a.fe, vcov=vcovHC (m214a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m214a.fe <- vcovHC (m214a.fe, method="white2", type="HC3", cluster="group"))

(sum.m214d.fe <- summary (m214d.fe <- plm (nLists.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m214d.fe <- coeftest (m214d.fe, vcov=vcovHC (m214d.fe, method="white2", type="HC3", cluster="group")))

# DV is ENPV
(sum.m224a.fe <- summary (m224a.fe <- plm (enpv.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m224a.fe <- coeftest (m224a.fe, vcov=vcovHC (m224a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m224a.fe <- vcovHC (m224a.fe, method="white2", type="HC3", cluster="group"))

(sum.m224d.fe <- summary (m224d.fe <- plm (enpv.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m224d.fe <- coeftest (m224d.fe, vcov=vcovHC (m224d.fe, method="white2", type="HC3", cluster="group")))

# DV is % vote for two largest parties
(sum.m234a.fe <- summary (m234a.fe <- plm (vote02sum.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m234a.fe <- coeftest (m234a.fe, vcov=vcovHC (m234a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m234a.fe <- vcovHC (m234a.fe, method="white2", type="HC3", cluster="group"))

(sum.m234d.fe <- summary (m234d.fe <- plm (vote02sum.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m234d.fe <- coeftest (m234d.fe, vcov=vcovHC (m234d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m234d.fe <- vcovHC (m234d.fe, method="white2", type="HC3", cluster="group"))



### Replicating the models reported in in Table A3 and Figure 4

## Panel (a): Main effect

# DV is # lists seats
(sum.m111a.fe <- summary (m111a.fe <- plm (nListSeats ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m111a.fe <- coeftest (m111a.fe, vcov=vcovHC (m111a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m111a.fe2 <- summary (m111a.fe2 <- plm (nListSeats ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))
(ci.m111a.fe2 <- cluster.wild.plm (m111a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m111d.fe <- summary (m111d.fe <- plm (nListSeats ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m111d.fe <- coeftest (m111d.fe, vcov=vcovHC (m111d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m111d.fe2 <- summary (m111d.fe2 <- plm (nListSeats ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m111d.fe2 <- cluster.wild.plm (m111d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

# DV is ENPS
(sum.m121a.fe <- summary (m121a.fe <- plm (enps ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m121a.fe <- coeftest (m121a.fe, vcov=vcovHC (m121a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m121a.fe2 <- summary (m121a.fe2 <- plm (enps ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))
(ci.m121a.fe2 <- cluster.wild.plm (m121a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m121d.fe <- summary (m121d.fe <- plm (enps ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m121d.fe <- coeftest (m121d.fe, vcov=vcovHC (m121d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m121d.fe2 <- summary (m121d.fe2 <- plm (enps ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m121d.fe2 <- cluster.wild.plm (m121d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

# DV is Gallagher index
(sum.m131a.fe <- summary (m131a.fe <- plm (dispGall ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m131a.fe <- coeftest (m131a.fe, vcov=vcovHC (m131a.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m131a.fe2 <- summary (m131a.fe2 <- plm (dispGall ~ largeMag + year, data=base, effect="individual", model="within", index=c ("district"))))
(ci.m131a.fe2 <- cluster.wild.plm (m131a.fe2, dat=base, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))

(sum.m131d.fe <- summary (m131d.fe <- plm (dispGall ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m131d.fe <- coeftest (m131d.fe, vcov=vcovHC (m131d.fe, method="white2", type="HC3", cluster="group")))

# wild-bootstrapped 95% CIs (NOTE: the bootstrapped SEs take a long time to compute)
(sum.m131d.fe2 <- summary (m131d.fe2 <- plm (dispGall ~ largeMag + year, data=base.sm, effect="individual", model="within", index=c ("district", "year"))))
(ci.m131d.fe2 <- cluster.wild.plm (m131d.fe2, dat=base.sm, cluster="group", ci.level = 0.95, boot.reps = B, report = TRUE, prog.bar = TRUE))


## Panel (b): Heterogeneous effects

# DV is # lists seats
(sum.m112a.fe <- summary (m112a.fe <- plm (nListSeats ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m112a.fe <- coeftest (m112a.fe, vcov=vcovHC (m112a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m112a.fe <- vcovHC (m112a.fe, method="white2", type="HC3", cluster="group"))

(sum.m112d.fe <- summary (m112d.fe <- plm (nListSeats ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m112d.fe <- coeftest (m112d.fe, vcov=vcovHC (m112d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m112d.fe <- vcovHC (m112d.fe, method="white2", type="HC3", cluster="group"))

# DV is ENPS
(sum.m122a.fe <- summary (m122a.fe <- plm (enps ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m122a.fe <- coeftest (m122a.fe, vcov=vcovHC (m122a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m122a.fe <- vcovHC (m122a.fe, method="white2", type="HC3", cluster="group"))

(sum.m122d.fe <- summary (m122d.fe <- plm (enps ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m122d.fe <- coeftest (m122d.fe, vcov=vcovHC (m122d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m122d.fe <- vcovHC (m122d.fe, method="white2", type="HC3", cluster="group"))

# DV is Gallagher index
(sum.m132a.fe <- summary (m132a.fe <- plm (dispGall ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m132a.fe <- coeftest (m132a.fe, vcov=vcovHC (m132a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m132a.fe <- vcovHC (m132a.fe, method="white2", type="HC3", cluster="group"))

(sum.m132d.fe <- summary (m132d.fe <- plm (dispGall ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m132d.fe <- coeftest (m132d.fe, vcov=vcovHC (m132d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m132d.fe <- vcovHC (m132d.fe, method="white2", type="HC3", cluster="group"))


## Panel (c): Rank-based DV

# DV is # lists seats
(sum.m113a.fe <- summary (m113a.fe <- plm (nListSeats.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m113a.fe <- coeftest (m113a.fe, vcov=vcovHC (m113a.fe, method="white2", type="HC3", cluster="group")))

(sum.m113d.fe <- summary (m113d.fe <- plm (nListSeats.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m113d.fe <- coeftest (m113d.fe, vcov=vcovHC (m113d.fe, method="white2", type="HC3", cluster="group")))

# DV is ENPS
(sum.m123a.fe <- summary (m123a.fe <- plm (enps.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m123a.fe <- coeftest (m123a.fe, vcov=vcovHC (m123a.fe, method="white2", type="HC3", cluster="group")))

(sum.m123d.fe <- summary (m123d.fe <- plm (enps.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m123d.fe <- coeftest (m123d.fe, vcov=vcovHC (m123d.fe, method="white2", type="HC3", cluster="group")))

# DV is Gallagher index
(sum.m133a.fe <- summary (m133a.fe <- plm (dispGall.rank ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m133a.fe <- coeftest (m133a.fe, vcov=vcovHC (m133a.fe, method="white2", type="HC3", cluster="group")))

(sum.m133d.fe <- summary (m133d.fe <- plm (dispGall.rank ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m133d.fe <- coeftest (m133d.fe, vcov=vcovHC (m133d.fe, method="white2", type="HC3", cluster="group")))


## Panel (d): Heterogeneous effects (rank)

# DV is # lists seats
(sum.m114a.fe <- summary (m114a.fe <- plm (nListSeats.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m114a.fe <- coeftest (m114a.fe, vcov=vcovHC (m114a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m114a.fe <- vcovHC (m114a.fe, method="white2", type="HC3", cluster="group"))

(sum.m114d.fe <- summary (m114d.fe <- plm (nListSeats.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m114d.fe <- coeftest (m114d.fe, vcov=vcovHC (m114d.fe, method="white2", type="HC3", cluster="group")))

# DV is ENPS
(sum.m124a.fe <- summary (m124a.fe <- plm (enps.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m124a.fe <- coeftest (m124a.fe, vcov=vcovHC (m124a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m124a.fe <- vcovHC (m124a.fe, method="white2", type="HC3", cluster="group"))

(sum.m124d.fe <- summary (m124d.fe <- plm (enps.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m124d.fe <- coeftest (m124d.fe, vcov=vcovHC (m124d.fe, method="white2", type="HC3", cluster="group")))

# DV is Gallagher index
(sum.m134a.fe <- summary (m134a.fe <- plm (dispGall.rank ~ largeMag + largeMag_vote3rd, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m134a.fe <- coeftest (m134a.fe, vcov=vcovHC (m134a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m134a.fe <- vcovHC (m134a.fe, method="white2", type="HC3", cluster="group"))

(sum.m134d.fe <- summary (m134d.fe <- plm (dispGall.rank ~ largeMag + largeMag_vote3rd, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m134d.fe <- coeftest (m134d.fe, vcov=vcovHC (m134d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m134d.fe <- vcovHC (m134d.fe, method="white2", type="HC3", cluster="group"))



### Replicating the models reported in Table A4 and Figure 5

## Panel (a): Main effect

# DV is # lists seats

# Mechanical effect: A -> B
(sum.m311a.fe <- summary (m311a.fe <- plm (nListSeats ~ largeAB, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m311a.fe <- coeftest (m311a.fe, vcov=vcovHC (m311a.fe, method="white2", type="HC3", cluster="group")))

(sum.m311d.fe <- summary (m311d.fe <- plm (nListSeats ~ largeAB, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m311d.fe <- coeftest (m311d.fe, vcov=vcovHC (m311d.fe, method="white2", type="HC3", cluster="group")))

# Psychological effect: B -> D
sum.m411a.fe <- summary (m411a.fe <- plm (nListSeats ~ largeBD, data=baseBD, effect="twoways", model="within", index=c ("district", "year2")))
(coef.m411a.fe <- coeftest (m411a.fe, vcov=vcovHC (m411a.fe, method="white2", type="HC3", cluster="group")))

(sum.m411d.fe <- summary (m411d.fe <- plm (nListSeats ~ largeBD, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m411d.fe <- coeftest (m411d.fe, vcov=vcovHC (m411d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (I): A -> C
(sum.m511a.fe <- summary (m511a.fe <- plm (nListSeats ~ largeAC, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m511a.fe <- coeftest (m511a.fe, vcov=vcovHC (m511a.fe, method="white2", type="HC3", cluster="group")))

(sum.m511d.fe <- summary (m511d.fe <- plm (nListSeats ~ largeAC, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m511d.fe <- coeftest (m511d.fe, vcov=vcovHC (m511d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (II): C -> D
(sum.m611a.fe <- summary (m611a.fe <- plm (nListSeats ~ largeCD, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m611a.fe <- coeftest (m611a.fe, vcov=vcovHC (m611a.fe, method="white2", type="HC3", cluster="group")))

(sum.m611d.fe <- summary (m611d.fe <- plm (nListSeats ~ largeCD, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m611d.fe <- coeftest (m611d.fe, vcov=vcovHC (m611d.fe, method="white2", type="HC3", cluster="group")))

# DV is ENPS

# Mechanical effect: A -> B
(sum.m321a.fe <- summary (m321a.fe <- plm (enps ~ largeAB, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m321a.fe <- coeftest (m321a.fe, vcov=vcovHC (m321a.fe, method="white2", type="HC3", cluster="group")))

(sum.m321d.fe <- summary (m321d.fe <- plm (enps ~ largeAB, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m321d.fe <- coeftest (m321d.fe, vcov=vcovHC (m321d.fe, method="white2", type="HC3", cluster="group")))

# Psychological effect: B -> D
(sum.m421a.fe <- summary (m421a.fe <- plm (enps ~ largeBD, data=baseBD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m421a.fe <- coeftest (m421a.fe, vcov=vcovHC (m421a.fe, method="white2", type="HC3", cluster="group")))

(sum.m421d.fe <- summary (m421d.fe <- plm (enps ~ largeBD, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m421d.fe <- coeftest (m421d.fe, vcov=vcovHC (m421d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (I): A -> C
(sum.m521a.fe <- summary (m521a.fe <- plm (enps ~ largeAC, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m521a.fe <- coeftest (m521a.fe, vcov=vcovHC (m521a.fe, method="white2", type="HC3", cluster="group")))

(sum.m521d.fe <- summary (m521d.fe <- plm (enps ~ largeAC, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m521d.fe <- coeftest (m521d.fe, vcov=vcovHC (m521d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (II): C -> D
(sum.m621a.fe <- summary (m621a.fe <- plm (enps ~ largeCD, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m621a.fe <- coeftest (m621a.fe, vcov=vcovHC (m621a.fe, method="white2", type="HC3", cluster="group")))

(sum.m621d.fe <- summary (m621d.fe <- plm (enps ~ largeCD, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m621d.fe <- coeftest (m621d.fe, vcov=vcovHC (m621d.fe, method="white2", type="HC3", cluster="group")))

# DV is Gallagher index

# Mechanical effect: A -> B
(sum.m331a.fe <- summary (m331a.fe <- plm (dispGall ~ largeAB, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m331a.fe <- coeftest (m331a.fe, vcov=vcovHC (m331a.fe, method="white2", type="HC3", cluster="group")))

(sum.m331d.fe <- summary (m331d.fe <- plm (dispGall ~ largeAB, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m331d.fe <- coeftest (m331d.fe, vcov=vcovHC (m331d.fe, method="white2", type="HC3", cluster="group")))

# Psychological effect: B -> D
(sum.m431a.fe <- summary (m431a.fe <- plm (dispGall ~ largeBD, data=baseBD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m431a.fe <- coeftest (m431a.fe, vcov=vcovHC (m431a.fe, method="white2", type="HC3", cluster="group")))

(sum.m431d.fe <- summary (m431d.fe <- plm (dispGall ~ largeBD, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m431d.fe <- coeftest (m431d.fe, vcov=vcovHC (m431d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (I): A -> C
(sum.m531a.fe <- summary (m531a.fe <- plm (dispGall ~ largeAC, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m531a.fe <- coeftest (m531a.fe, vcov=vcovHC (m531a.fe, method="white2", type="HC3", cluster="group")))

(sum.m531d.fe <- summary (m531d.fe <- plm (dispGall ~ largeAC, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m531d.fe <- coeftest (m531d.fe, vcov=vcovHC (m531d.fe, method="white2", type="HC3", cluster="group")))

# Component psychological effect (II): C -> D
(sum.m631a.fe <- summary (m631a.fe <- plm (dispGall ~ largeCD, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m631a.fe <- coeftest (m631a.fe, vcov=vcovHC (m631a.fe, method="white2", type="HC3", cluster="group")))

(sum.m631d.fe <- summary (m631d.fe <- plm (dispGall ~ largeCD, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m631d.fe <- coeftest (m631d.fe, vcov=vcovHC (m631d.fe, method="white2", type="HC3", cluster="group")))


## Panel (b): Heterogeneous effect

# DV is # lists seats

# Mechanical effect: A -> B
(sum.m312a.fe <- summary (m312a.fe <- plm (nListSeats ~ largeAB + largeAB_vote3rd, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m312a.fe <- coeftest (m312a.fe, vcov=vcovHC (m312a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m312a.fe <- vcovHC (m312a.fe, method="white2", type="HC3", cluster="group"))

(sum.m312d.fe <- summary (m312d.fe <- plm (nListSeats ~ largeAB + largeAB_vote3rd, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m312d.fe <- coeftest (m312d.fe, vcov=vcovHC (m312d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m312d.fe <- vcovHC (m312d.fe, method="white2", type="HC3", cluster="group"))

# Psychological effect: B -> D
(sum.m412a.fe <- summary (m412a.fe <- plm (nListSeats ~ largeBD + largeBD_vote3rd, data=baseBD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m412a.fe <- coeftest (m412a.fe, vcov=vcovHC (m412a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m412a.fe <- vcovHC (m412a.fe, method="white2", type="HC3", cluster="group"))

(sum.m412d.fe <- summary (m412d.fe <- plm (nListSeats ~ largeBD + largeBD_vote3rd, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m412d.fe <- coeftest (m412d.fe, vcov=vcovHC (m412d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m412d.fe <- vcovHC (m412d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (I): A -> C
(sum.m512a.fe <- summary (m512a.fe <- plm (nListSeats ~ largeAC + largeAC_vote3rd, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m512a.fe <- coeftest (m512a.fe, vcov=vcovHC (m512a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m512a.fe <- vcovHC (m512a.fe, method="white2", type="HC3", cluster="group"))

(sum.m512d.fe <- summary (m512d.fe <- plm (nListSeats ~ largeAC + largeAC_vote3rd, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m512d.fe <- coeftest (m512d.fe, vcov=vcovHC (m512d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m512d.fe <- vcovHC (m512d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (II): C -> D
(sum.m612a.fe <- summary (m612a.fe <- plm (nListSeats ~ largeCD + largeCD_vote3rd, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m612a.fe <- coeftest (m612a.fe, vcov=vcovHC (m612a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m612a.fe <- vcovHC (m612a.fe, method="white2", type="HC3", cluster="group"))

(sum.m612d.fe <- summary (m612d.fe <- plm (nListSeats ~ largeCD + largeCD_vote3rd, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m612d.fe <- coeftest (m612d.fe, vcov=vcovHC (m612d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m612d.fe <- vcovHC (m612d.fe, method="white2", type="HC3", cluster="group"))

# DV is ENPS

# Mechanical effect: A -> B
(sum.m322a.fe <- summary (m322a.fe <- plm (enps ~ largeAB + largeAB_vote3rd, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m322a.fe <- coeftest (m322a.fe, vcov=vcovHC (m322a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m322a.fe <- vcovHC (m322a.fe, method="white2", type="HC3", cluster="group"))

(sum.m322d.fe <- summary (m322d.fe <- plm (enps ~ largeAB + largeAB_vote3rd, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m322d.fe <- coeftest (m322d.fe, vcov=vcovHC (m322d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m322d.fe <- vcovHC (m322d.fe, method="white2", type="HC3", cluster="group"))

# Psychological effect: B -> D
(sum.m422a.fe <- summary (m422a.fe <- plm (enps ~ largeBD + largeBD_vote3rd, data=baseBD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m422a.fe <- coeftest (m422a.fe, vcov=vcovHC (m422a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m422a.fe <- vcovHC (m422a.fe, method="white2", type="HC3", cluster="group"))

(sum.m422d.fe <- summary (m422d.fe <- plm (enps ~ largeBD + largeBD_vote3rd, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m422d.fe <- coeftest (m422d.fe, vcov=vcovHC (m422d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m422d.fe <- vcovHC (m422d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (I): A -> C
(sum.m522a.fe <- summary (m522a.fe <- plm (enps ~ largeAC + largeAC_vote3rd, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m522a.fe <- coeftest (m522a.fe, vcov=vcovHC (m522a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m522a.fe <- vcovHC (m522a.fe, method="white2", type="HC3", cluster="group"))

(sum.m522d.fe <- summary (m522d.fe <- plm (enps ~ largeAC + largeAC_vote3rd, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m522d.fe <- coeftest (m522d.fe, vcov=vcovHC (m522d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m522d.fe <- vcovHC (m522d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (II): C -> D
(sum.m622a.fe <- summary (m622a.fe <- plm (enps ~ largeCD + largeCD_vote3rd, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m622a.fe <- coeftest (m622a.fe, vcov=vcovHC (m622a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m622a.fe <- vcovHC (m622a.fe, method="white2", type="HC3", cluster="group"))

(sum.m622d.fe <- summary (m622d.fe <- plm (enps ~ largeCD + largeCD_vote3rd, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m622d.fe <- coeftest (m622d.fe, vcov=vcovHC (m622d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m622d.fe <- vcovHC (m622d.fe, method="white2", type="HC3", cluster="group"))

# DV is Gallagher index

# Mechanical effect: A -> B
(sum.m332a.fe <- summary (m332a.fe <- plm (dispGall ~ largeAB + largeAB_vote3rd, data=baseAB, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m332a.fe <- coeftest (m332a.fe, vcov=vcovHC (m332a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m332a.fe <- vcovHC (m332a.fe, method="white2", type="HC3", cluster="group"))

(sum.m332d.fe <- summary (m332d.fe <- plm (dispGall ~ largeAB + largeAB_vote3rd, data=baseAB.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m332d.fe <- coeftest (m332d.fe, vcov=vcovHC (m332d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m332d.fe <- vcovHC (m332d.fe, method="white2", type="HC3", cluster="group"))

# Psychological effect: B -> D
(sum.m432a.fe <- summary (m432a.fe <- plm (dispGall ~ largeBD + largeBD_vote3rd, data=baseBD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m432a.fe <- coeftest (m432a.fe, vcov=vcovHC (m432a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m432a.fe <- vcovHC (m432a.fe, method="white2", type="HC3", cluster="group"))

(sum.m432d.fe <- summary (m432d.fe <- plm (dispGall ~ largeBD + largeBD_vote3rd, data=baseBD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m432d.fe <- coeftest (m432d.fe, vcov=vcovHC (m432d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m432d.fe <- vcovHC (m432d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (I): A -> C
(sum.m532a.fe <- summary (m532a.fe <- plm (dispGall ~ largeAC + largeAC_vote3rd, data=baseAC, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m532a.fe <- coeftest (m532a.fe, vcov=vcovHC (m532a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m532a.fe <- vcovHC (m532a.fe, method="white2", type="HC3", cluster="group"))

(sum.m532d.fe <- summary (m532d.fe <- plm (dispGall ~ largeAC + largeAC_vote3rd, data=baseAC.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m532d.fe <- coeftest (m532d.fe, vcov=vcovHC (m532d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m532d.fe <- vcovHC (m532d.fe, method="white2", type="HC3", cluster="group"))

# Component psychological effect (II): C -> D
(sum.m632a.fe <- summary (m632a.fe <- plm (dispGall ~ largeCD + largeCD_vote3rd, data=baseCD, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m632a.fe <- coeftest (m632a.fe, vcov=vcovHC (m632a.fe, method="white2", type="HC3", cluster="group")))
(vcov.m632a.fe <- vcovHC (m632a.fe, method="white2", type="HC3", cluster="group"))

(sum.m632d.fe <- summary (m632d.fe <- plm (dispGall ~ largeCD + largeCD_vote3rd, data=baseCD.sm, effect="twoways", model="within", index=c ("district", "year2"))))
(coef.m632d.fe <- coeftest (m632d.fe, vcov=vcovHC (m632d.fe, method="white2", type="HC3", cluster="group")))
(vcov.m632d.fe <- vcovHC (m632d.fe, method="white2", type="HC3", cluster="group"))



# Bootstrapping the difference between coefficients: [C -> D] - [A -> B]

# creating two "short" datasets
baseFull2 <- baseFull[baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),]
baseFull2$district <- factor (baseFull2$district)

baseFull2.sm <- baseFull[baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,]
baseFull2.sm$district <- factor (baseFull2.sm$district)

baseFull2$districtB <- NA; baseFull2.sm$districtB <- NA

# creating the datasets to "fill-in" later
coef.full <- coef.sm <- matrix (NA, nrow=B, ncol=18)

# doing the bootstrap
set.seed (193974) ## for reproducibility

# WARNING! The bootstrapped is commented out because it took approx. 23 minutes to run
# system.time (
#   for (b in 1:B){
#     
#     ## sampling provinces
#     G <- 19; G.sm <- 10
#     samp.full <- sample (unique (baseFull2$district), size=G, replace=T)
#     samp.sm <- sample (unique (baseFull2.sm$district), size=G.sm, replace=T)
#     
#     ## creating the new datasets
#     baseFull.tmp <- NULL
#     baseFull.sm.tmp <- NULL
#     
#     for (i in 1:G){
#       baseFull2[baseFull2$district==samp.full[i],]$districtB <- paste ("district", i, sep="_")
#       baseFull.tmp <- rbind.fill (baseFull.tmp, baseFull2[baseFull2$district==samp.full[i],])
#     }
#     baseFull.tmp$districtB <- factor (baseFull.tmp$districtB)
#     baseFull.tmp$ffCode <- factor (baseFull.tmp$ffCode)
#     
#     for (i in 1:G.sm){
#       baseFull2.sm[baseFull2.sm$district==samp.sm[i],]$districtB <- paste ("district", i, sep="_")
#       baseFull.sm.tmp <- rbind.fill (baseFull.sm.tmp, baseFull2.sm[baseFull2.sm$district==samp.sm[i],])
#     }
#     baseFull.sm.tmp$districtB <- factor (baseFull.sm.tmp$districtB)
#     baseFull.sm.tmp$ffCode <- factor (baseFull.sm.tmp$ffCode)
#     
#     # .plm versions
#     baseAB.tmp <- pdata.frame (baseFull.tmp[!is.na (baseFull.tmp$ffCode) & (baseFull.tmp$ffCode=="A" | baseFull.tmp$ffCode=="B"),], index=c ("districtB", "year2"))
#     baseCD.tmp <- pdata.frame (baseFull.tmp[!is.na (baseFull.tmp$ffCode) & (baseFull.tmp$ffCode=="C" | baseFull.tmp$ffCode=="D"),], index=c ("districtB", "year2"))
#     
#     baseAB.sm.tmp <- pdata.frame (baseFull.sm.tmp[!is.na (baseFull.sm.tmp$ffCode) & (baseFull.sm.tmp$ffCode=="A" | baseFull.sm.tmp$ffCode=="B"),], index=c ("districtB", "year2"))
#     baseCD.sm.tmp <- pdata.frame (baseFull.sm.tmp[!is.na (baseFull.sm.tmp$ffCode) & (baseFull.sm.tmp$ffCode=="C" | baseFull.sm.tmp$ffCode=="D"),], index=c ("districtB", "year2"))
#     
#     ## fitting the models
#     
#     # (a) full sample
#     
#     # number of parties getting seats
#     coef.full[b,1] <- plm (nListSeats ~ largeCD, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,2] <- plm (nListSeats ~ largeAB, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,3:4] <- plm (nListSeats ~ largeCD + largeCD_vote3rd, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,5:6] <- plm (nListSeats ~ largeAB + largeAB_vote3rd, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     
#     # ENPS
#     coef.full[b,7] <- plm (enps ~ largeCD, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,8] <- plm (enps ~ largeAB, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,9:10] <- plm (enps ~ largeCD + largeCD_vote3rd, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,11:12] <- plm (enps ~ largeAB + largeAB_vote3rd, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     
#     # Gallagher Index  
#     coef.full[b,13] <- plm (dispGall ~ largeCD, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,14] <- plm (dispGall ~ largeAB, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,15:16] <- plm (dispGall ~ largeCD + largeCD_vote3rd, data=baseCD.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.full[b,17:18] <- plm (dispGall ~ largeAB + largeAB_vote3rd, data=baseAB.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     
#     
#     # (b) small-province sample
#     
#     # number of parties getting seats
#     coef.sm[b,1] <- plm (nListSeats ~ largeCD, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,2] <- plm (nListSeats ~ largeAB, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,3:4] <- plm (nListSeats ~ largeCD + largeCD_vote3rd, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,5:6] <- plm (nListSeats ~ largeAB + largeAB_vote3rd, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     
#     # ENPS
#     coef.sm[b,7] <- plm (enps ~ largeCD, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,8] <- plm (enps ~ largeAB, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,9:10] <- plm (enps ~ largeCD + largeCD_vote3rd, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,11:12] <- plm (enps ~ largeAB + largeAB_vote3rd, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     
#     # Gallagher Index  
#     coef.sm[b,13] <- plm (dispGall ~ largeCD, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,14] <- plm (dispGall ~ largeAB, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,15:16] <- plm (dispGall ~ largeCD + largeCD_vote3rd, data=baseCD.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#     coef.sm[b,17:18] <- plm (dispGall ~ largeAB + largeAB_vote3rd, data=baseAB.sm.tmp, effect="twoways", model="within", index=c ("districtB", "year2"))$coef
#   }
# ) ## 23 min (1358 seconds)



### Replicating the models reported in in Table A5 and Figure 6

# DV is revenues per capita (log)
(sum.m811a.fe <- summary (m811a.fe <- plm (log (revenues.head) ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m811a.fe <- coeftest (m811a.fe, vcov=vcovHC (m811a.fe, method="white2", type="HC3", cluster="group")))

(sum.m811b.fe <- summary (m811b.fe <- plm (log (revenues.head) ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m811b.fe <- coeftest (m811b.fe, vcov=vcovHC (m811b.fe, method="white2", type="HC3", cluster="group")))

# DV is % own revenues
(sum.m821a.fe <- summary (m821a.fe <- plm (rev.own ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m821a.fe <- coeftest (m821a.fe, vcov=vcovHC (m821a.fe, method="white2", type="HC3", cluster="group")))

(sum.m821b.fe <- summary (m821b.fe <- plm (rev.own ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m821b.fe <- coeftest (m821b.fe, vcov=vcovHC (m821b.fe, method="white2", type="HC3", cluster="group")))

# DV is % royalties
(sum.m831a.fe <- summary (m831a.fe <- plm (rev.royalties ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m831a.fe <- coeftest (m831a.fe, vcov=vcovHC (m831a.fe, method="white2", type="HC3", cluster="group")))

(sum.m831b.fe <- summary (m831b.fe <- plm (rev.royalties ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m831b.fe <- coeftest (m831b.fe, vcov=vcovHC (m831b.fe, method="white2", type="HC3", cluster="group")))

# DV is % automatic transfers
(sum.m841a.fe <- summary (m841a.fe <- plm (rev.grant.auto ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m841a.fe <- coeftest (m841a.fe, vcov=vcovHC (m841a.fe, method="white2", type="HC3", cluster="group")))

(sum.m841b.fe <- summary (m841b.fe <- plm (rev.grant.auto ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m841b.fe <- coeftest (m841b.fe, vcov=vcovHC (m841b.fe, method="white2", type="HC3", cluster="group")))

# DV is % discretionary transfers
(sum.m851a.fe <- summary (m851a.fe <- plm (rev.grant.discret ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m851a.fe <- coeftest (m851a.fe, vcov=vcovHC (m851a.fe, method="white2", type="HC3", cluster="group")))

(sum.m851b.fe <- summary (m851b.fe <- plm (rev.grant.discret ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m851b.fe <- coeftest (m851b.fe, vcov=vcovHC (m851b.fe, method="white2", type="HC3", cluster="group")))

# DV is public employees per 1,000 inhabitants
(sum.m861a.fe <- summary (m861a.fe <- plm (employees.th ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m861a.fe <- coeftest (m861a.fe, vcov=vcovHC (m861a.fe, method="white2", type="HC3", cluster="group")))

(sum.m861b.fe <- summary (m861b.fe <- plm (employees.th ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m861b.fe <- coeftest (m861b.fe, vcov=vcovHC (m861b.fe, method="white2", type="HC3", cluster="group")))

# DV is unemployment rate
(sum.m871a.fe <- summary (m871a.fe <- plm (unemployment ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m871a.fe <- coeftest (m871a.fe, vcov=vcovHC (m871a.fe, method="white2", type="HC3", cluster="group")))

(sum.m871b.fe <- summary (m871b.fe <- plm (unemployment ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m871b.fe <- coeftest (m871b.fe, vcov=vcovHC (m871b.fe, method="white2", type="HC3", cluster="group")))

# DV is infant mortality rate
(sum.m881a.fe <- summary (m881a.fe <- plm (imr ~ largeMag, data=base, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m881a.fe <- coeftest (m881a.fe, vcov=vcovHC (m881a.fe, method="white2", type="HC3", cluster="group")))

(sum.m881b.fe <- summary (m881b.fe <- plm (imr ~ largeMag, data=base.sm, effect="twoways", model="within", index=c ("district", "year"))))
(coef.m881b.fe <- coeftest (m881b.fe, vcov=vcovHC (m881b.fe, method="white2", type="HC3", cluster="group")))




########## BUILDING THE TABLES

## Setting values that we will keep constant when presenting the results
(digits <- 3) ## number of decimals to use
v.high <- 15; v.low <- 5 ## high and low values for heterogeneous effects
(tval <- qt (.975, df=18)); (tval.sm <- qt (.975, df=9))  ## for buildings 95% CIs


### Replicating Table 2: Descriptive statistics
(tab.descriptive01 <- cbind (
  
  rep (NA, 8) ##empty column
  
  ## Full sample
  , with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], rbind (
    
    # Explanatory variables
    c (mean (mag, na.rm=TRUE), sd (mag, na.rm=TRUE), min (mag, na.rm=TRUE), max (mag, na.rm=TRUE))
    , c (mean (vote3rd, na.rm=TRUE), sd (vote3rd, na.rm=TRUE), min (vote3rd, na.rm=TRUE), max (vote3rd, na.rm=TRUE))
    
    # DV: Electoral coordination
    , c (mean (nLists, na.rm=TRUE), sd (nLists, na.rm=TRUE), min (nLists, na.rm=TRUE), max (nLists, na.rm=TRUE))
    , c (mean (enpv, na.rm=TRUE), sd (enpv, na.rm=TRUE), min (enpv, na.rm=TRUE), max (enpv, na.rm=TRUE))
    , c (mean (vote02sum, na.rm=TRUE), sd (vote02sum, na.rm=TRUE), min (vote02sum, na.rm=TRUE), max (vote02sum, na.rm=TRUE))
    
    # DV: Aggregate outcomes
    , c (mean (nListSeats, na.rm=TRUE), sd (nListSeats, na.rm=TRUE), min (nListSeats, na.rm=TRUE), max (nListSeats, na.rm=TRUE))
    , c (mean (enps, na.rm=TRUE), sd (enps, na.rm=TRUE), min (enps, na.rm=TRUE), max (enps, na.rm=TRUE))
    , c (mean (dispGall, na.rm=TRUE), sd (dispGall, na.rm=TRUE), min (dispGall, na.rm=TRUE), max (dispGall, na.rm=TRUE))
  ))
  
  , rep (NA, 8)
  
  ## Small provinces (Delegation <= 5)
  , with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], rbind (
    
    # Explanatory variables
    c (mean (mag, na.rm=TRUE), sd (mag, na.rm=TRUE), min (mag, na.rm=TRUE), max (mag, na.rm=TRUE))
    , c (mean (vote3rd, na.rm=TRUE), sd (vote3rd, na.rm=TRUE), min (vote3rd, na.rm=TRUE), max (vote3rd, na.rm=TRUE))
    
    # DV: Electoral coordination
    , c (mean (nLists, na.rm=TRUE), sd (nLists, na.rm=TRUE), min (nLists, na.rm=TRUE), max (nLists, na.rm=TRUE))
    , c (mean (enpv, na.rm=TRUE), sd (enpv, na.rm=TRUE), min (enpv, na.rm=TRUE), max (enpv, na.rm=TRUE))
    , c (mean (vote02sum, na.rm=TRUE), sd (vote02sum, na.rm=TRUE), min (vote02sum, na.rm=TRUE), max (vote02sum, na.rm=TRUE))
    
    # DV: Aggregate outcomes
    , c (mean (nListSeats, na.rm=TRUE), sd (nListSeats, na.rm=TRUE), min (nListSeats, na.rm=TRUE), max (nListSeats, na.rm=TRUE))
    , c (mean (enps, na.rm=TRUE), sd (enps, na.rm=TRUE), min (enps, na.rm=TRUE), max (enps, na.rm=TRUE))
    , c (mean (dispGall, na.rm=TRUE), sd (dispGall, na.rm=TRUE), min (dispGall, na.rm=TRUE), max (dispGall, na.rm=TRUE))
  ))
))

## within-province standard deviations
(sd.all <- c (
  mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (mag, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], sd (vote3rd, na.rm=TRUE)))
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (nLists, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (enpv, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (vote02sum, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (nListSeats, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (enps, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),], by (dispGall, district, sd, na.rm=TRUE))[], na.rm=T)
))
(sd.sm <- c (
  mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (mag, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], sd (vote3rd, na.rm=TRUE)))
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (nLists, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (enpv, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (vote02sum, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (nListSeats, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (enps, district, sd, na.rm=TRUE))[], na.rm=T)
  , mean (with (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,], by (dispGall, district, sd, na.rm=TRUE))[], na.rm=T)
))
tab.descriptive01[,3] <- sd.all
tab.descriptive01[,8] <- sd.sm

rows <- c (
  "\\emph{magnitude}", "\\emph{vote third party}" ## explanatory variables
  , "\\emph{\\# lists running}", "\\emph{ENPV}", "\\emph{vote first two}"  ## DV (2)
  , "\\emph{\\# lists seats}", "\\emph{ENPS}", "\\emph{Gallagher index}"  ## DV (1)
  )

(n.full <- nrow (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989),]))
(n.sm <- nrow (baseFull[baseFull$status=="actual" & baseFull$year!=1983 & baseFull$year!=1994 & baseFull$odd==1 & !(baseFull$district=="TDF" & baseFull$year < 1989) & baseFull$delegation <= 5,]))

Header1 <- paste ("\\toprule \\toprule & & \\multicolumn{4}{c}{\\textbf{Full sample}} & & \\multicolumn{4}{c}{\\textbf{Small provinces}} \\\\ \n")
Header2 <- paste ("& & \\multicolumn{4}{c}{($\\text{\\# provinces: $19$; }n = ", sprintf ("%.0f", round (n.full,0)), "$)} & & \\multicolumn{4}{c}{($\\text{\\# provinces: $10$; }n = ", sprintf ("%.0f", round (n.sm,0)), "$)} \\\\ \\cmidrule{3-6} \\cmidrule{8-11} \n")
Header3 <- paste ("\\multicolumn{11}{l}{\\textbf{(a) Explanatory}} \\\\ \n")
Header4 <- paste ("\\multicolumn{2}{l}{~~~~~\\textbf{variables}} & \\multicolumn{1}{c}{mean} & \\multicolumn{1}{c}{sd.$^*$} & \\multicolumn{1}{c}{min.} & \\multicolumn{1}{c}{max.} & & \\multicolumn{1}{c}{mean} & \\multicolumn{1}{c}{sd.$^*$} & \\multicolumn{1}{c}{min.} & \\multicolumn{1}{c}{max.} \\\\ \\midrule \n")
Header5 <- paste ("[2.5ex] \\multicolumn{11}{l}{\\textbf{(b) Dependent variables (1): Electoral coordination}} \\\\ \\midrule \n")
Header6 <- paste ("[2.5ex] \\multicolumn{11}{l}{\\textbf{(c) Dependent variables (2): Seat distribution}} \\\\ \\midrule \n")

addtorow <- list ()
addtorow$pos <- list ()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 0
addtorow$pos[[5]] <- 2
addtorow$pos[[6]] <- 5
addtorow$command <- c (Header1, Header2, Header3, Header4
                       , Header5, Header6)
print (xtable ( cbind (rows, as.data.frame (tab.descriptive01))
                , align=c ("l","l","r","r","r","r","r","c","r","r","r","r")
                , digits=digits-1
                , caption="\\emph{Descriptive statistics}"
                , label="T:descriptive")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="small"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )


### Replicating Table A2: Electoral coordination
(tab02coef <- cbind (
  c (coef.m211a.fe[1,1], coef.m212a.fe[1:2,1], coef.m213a.fe[1,1], coef.m214a.fe[1:2,1])
  , c (coef.m211d.fe[1,1], coef.m212d.fe[1:2,1], coef.m213d.fe[1,1], coef.m214d.fe[1:2,1])
  , c (coef.m221a.fe[1,1], coef.m222a.fe[1:2,1], coef.m223a.fe[1,1], coef.m224a.fe[1:2,1])
  , c (coef.m221d.fe[1,1], coef.m222d.fe[1:2,1], coef.m223d.fe[1,1], coef.m224d.fe[1:2,1])
  , c (coef.m231a.fe[1,1], coef.m232a.fe[1:2,1], coef.m233a.fe[1,1], coef.m234a.fe[1:2,1])
  , c (coef.m231d.fe[1,1], coef.m232d.fe[1:2,1], coef.m233d.fe[1,1], coef.m234d.fe[1:2,1])
))
(tab02se <- cbind (
  c (coef.m211a.fe[1,2], coef.m212a.fe[1:2,2], coef.m213a.fe[1,2], coef.m214a.fe[1:2,2])
  , c (coef.m211d.fe[1,2], coef.m212d.fe[1:2,2], coef.m213d.fe[1,2], coef.m214d.fe[1:2,2])
  , c (coef.m221a.fe[1,2], coef.m222a.fe[1:2,2], coef.m223a.fe[1,2], coef.m224a.fe[1:2,2])
  , c (coef.m221d.fe[1,2], coef.m222d.fe[1:2,2], coef.m223d.fe[1,2], coef.m224d.fe[1:2,2])
  , c (coef.m231a.fe[1,2], coef.m232a.fe[1:2,2], coef.m233a.fe[1,2], coef.m234a.fe[1:2,2])
  , c (coef.m231d.fe[1,2], coef.m232d.fe[1:2,2], coef.m233d.fe[1,2], coef.m234d.fe[1:2,2])
))

# estimates with 95% CIs below
(tab02ci.full <- paste ("["
                        , sprintf ("%.2f", round (c (tab02coef[,c (1,3,5)]) + c (tab02se[,c (1,3,5)])*(-1)*tval, digits))
                        , ":"
                        , sprintf ("%.2f", round (c (tab02coef[,c (1,3,5)]) + c (tab02se[,c (1,3,5)])*(1)*tval, digits))
                        , "]", sep=""))
(tab02ci.sm <- paste ("["
                      , sprintf ("%.2f", round (c (tab02coef[,c (2,4,6)]) + c (tab02se[,c (2,4,6)])*(-1)*tval, digits))
                      , ":"
                      , sprintf ("%.2f", round (c (tab02coef[,c (2,4,6)]) + c (tab02se[,c (2,4,6)])*(1)*tval, digits))
                      , "]", sep=""))
(tab02ci.wild <- paste ("["
                        , sprintf ("%.2f", round (c (ci.m211a.fe2$ci[1,1], ci.m211d.fe2$ci[1,1], ci.m221a.fe2$ci[1,1], ci.m221d.fe2$ci[1,1], ci.m231a.fe2$ci[1,1], ci.m231d.fe2$ci[1,1]), digits))
                        , ":"
                        , sprintf ("%.2f", round (c (ci.m211a.fe2$ci[1,2], ci.m211d.fe2$ci[1,2], ci.m221a.fe2$ci[1,2], ci.m221d.fe2$ci[1,2], ci.m231a.fe2$ci[1,2], ci.m231d.fe2$ci[1,2]), digits))
                        , "]", sep=""))

# creating the table
odd <- seq (1, 12, by=2); even <- seq (2,12, by=2)
(tab02all <- rbind (
  matrix (sprintf ("%.2f", round (tab02coef, digits)), nrow=6, byrow=F)
  , matrix (c (tab02ci.full, tab02ci.sm)[c (1:6, 19:24, 7:12, 25:30, 13:18, 31:36)], nrow=6, byrow=F)
)[sort (rep (1:6, 2)) + rep (c (0, 6), 6),])
(tab02all <- rbind (tab02all[c (1:2),], tab02ci.wild, tab02all[c (3:nrow (tab02all)),]))

rows <- c (
  "\\emph{magnitude}", "", ""
  , "\\emph{magnitude}", "", "\\emph{magnitude}", "~~~~~$\\times$ \\emph{vote third party}"
  , "\\emph{magnitude}", ""
  , "\\emph{magnitude}", "", "\\emph{magnitude}", "~~~~~$\\times$ \\emph{vote third party}")

(n.obs <- cbind (unlist (pdim (m211a.fe)$nT), unlist (pdim (m211d.fe)$nT), unlist (pdim (m221a.fe)$nT), unlist (pdim (m221d.fe)$nT), unlist (pdim (m231a.fe)$nT), unlist (pdim (m231d.fe)$nT))[c (3,1,2),])

Header1 <- paste ("\\toprule \\toprule & \\multicolumn{2}{c}{\\textbf{\\emph{\\# lists running}}} & \\multicolumn{2}{c}{\\textbf{\\emph{ENPV}}} & \\multicolumn{2}{c}{\\textbf{\\emph{vote first two}}} \\\\ \n")
Header2 <- paste ("[0.5ex] \\multicolumn{1}{l}{} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small}& \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} \\\\ \n")
Header3 <- paste ("\\multicolumn{1}{l}{\\textbf{(a) Main effect}} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(b) Heterogeneous effect}} \\\\ \\midrule \n")
Header5 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(c) Main effect (rank)}} \\\\ \\midrule \n")
Header6 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(b) Heterogeneous effect (rank)}} \\\\ \\midrule \n")

add.nobs <- paste ("[2.5ex] \\midrule num. obs & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,6],0)), "} \\\\ \n")
add.provs <- paste ("provinces & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,6],0)), "} \\\\ \n")
add.time <- paste ("elections & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,6],0)), "} \\\\ \n")

addtorow <- list ()
addtorow$pos <- list ()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 3
addtorow$pos[[5]] <- 7
addtorow$pos[[6]] <- 9
addtorow$pos[[7]] <- 13
addtorow$pos[[8]] <- 13
addtorow$pos[[9]] <- 13
addtorow$command <- c (Header1, Header2, Header3
                       , Header4, Header5, Header6
                       , add.nobs, add.provs, add.time)
print (xtable ( cbind (rows, tab02all)
                , align=c ("l","l","c","c","c","c","c","c")
                , digits=digits
                , caption="\\emph{District magnitude and electoral coordination in Argentina, 1985-2015}"
                , label="T:results02")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="small"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )



### Replicating Table A3: Aggregate outcomes
(tab01coef <- cbind (
  c (coef.m111a.fe[1,1], coef.m112a.fe[1:2,1], coef.m113a.fe[1,1], coef.m114a.fe[1:2,1])
  , c (coef.m111d.fe[1,1], coef.m112d.fe[1:2,1], coef.m113d.fe[1,1], coef.m114d.fe[1:2,1])
  , c (coef.m121a.fe[1,1], coef.m122a.fe[1:2,1], coef.m123a.fe[1,1], coef.m124a.fe[1:2,1])
  , c (coef.m121d.fe[1,1], coef.m122d.fe[1:2,1], coef.m123d.fe[1,1], coef.m124d.fe[1:2,1])
  , c (coef.m131a.fe[1,1], coef.m132a.fe[1:2,1], coef.m133a.fe[1,1], coef.m134a.fe[1:2,1])
  , c (coef.m131d.fe[1,1], coef.m132d.fe[1:2,1], coef.m133d.fe[1,1], coef.m134d.fe[1:2,1])
))
(tab01se <- cbind (
  c (coef.m111a.fe[1,2], coef.m112a.fe[1:2,2], coef.m113a.fe[1,2], coef.m114a.fe[1:2,2])
  , c (coef.m111d.fe[1,2], coef.m112d.fe[1:2,2], coef.m113d.fe[1,2], coef.m114d.fe[1:2,2])
  , c (coef.m121a.fe[1,2], coef.m122a.fe[1:2,2], coef.m123a.fe[1,2], coef.m124a.fe[1:2,2])
  , c (coef.m121d.fe[1,2], coef.m122d.fe[1:2,2], coef.m123d.fe[1,2], coef.m124d.fe[1:2,2])
  , c (coef.m131a.fe[1,2], coef.m132a.fe[1:2,2], coef.m133a.fe[1,2], coef.m134a.fe[1:2,2])
  , c (coef.m131d.fe[1,2], coef.m132d.fe[1:2,2], coef.m133d.fe[1,2], coef.m134d.fe[1:2,2])
))

# estimates with 95% CIs below
(tab01ci.full <- paste ("["
                        , sprintf ("%.2f", round (c (tab01coef[,c (1,3,5)]) + c (tab01se[,c (1,3,5)])*(-1)*tval, digits))
                        , ":"
                        , sprintf ("%.2f", round (c (tab01coef[,c (1,3,5)]) + c (tab01se[,c (1,3,5)])*(1)*tval, digits))
                        , "]", sep=""))
(tab01ci.sm <- paste ("["
                      , sprintf ("%.2f", round (c (tab01coef[,c (2,4,6)]) + c (tab01se[,c (2,4,6)])*(-1)*tval, digits))
                      , ":"
                      , sprintf ("%.2f", round (c (tab01coef[,c (2,4,6)]) + c (tab01se[,c (2,4,6)])*(1)*tval, digits))
                      , "]", sep=""))
(tab01ci.wild <- paste ("["
                        , sprintf ("%.2f", round (c (ci.m111a.fe2$ci[1,1], ci.m111d.fe2$ci[1,1], ci.m121a.fe2$ci[1,1], ci.m121d.fe2$ci[1,1], ci.m131a.fe2$ci[1,1], ci.m131d.fe2$ci[1,1]), digits))
                        , ":"
                        , sprintf ("%.2f", round (c (ci.m111a.fe2$ci[1,2], ci.m111d.fe2$ci[1,2], ci.m121a.fe2$ci[1,2], ci.m121d.fe2$ci[1,2], ci.m131a.fe2$ci[1,2], ci.m131d.fe2$ci[1,2]), digits))
                        , "]", sep=""))

# creating the table
odd <- seq (1, 12, by=2); even <- seq (2,12, by=2)
(tab01all <- rbind (
  matrix (sprintf ("%.2f", round (tab01coef, digits)), nrow=6, byrow=F)
  , matrix (c (tab01ci.full, tab01ci.sm)[c (1:6, 19:24, 7:12, 25:30, 13:18, 31:36)], nrow=6, byrow=F)
)[sort (rep (1:6, 2)) + rep (c (0, 6), 6),])
(tab01all <- rbind (tab01all[c (1:2),], tab01ci.wild, tab01all[c (3:nrow (tab01all)),]))

rows <- c (
  "\\emph{magnitude}", "", ""
  , "\\emph{magnitude}", "", "\\emph{magnitude}", "~~~~~$\\times$ \\emph{vote third party}"
  , "\\emph{magnitude}", ""
  , "\\emph{magnitude}", "", "\\emph{magnitude}", "~~~~~$\\times$ \\emph{vote third party}")

(n.obs <- cbind (unlist (pdim (m111a.fe)$nT), unlist (pdim (m111d.fe)$nT), unlist (pdim (m121a.fe)$nT), unlist (pdim (m121d.fe)$nT), unlist (pdim (m131a.fe)$nT), unlist (pdim (m131d.fe)$nT))[c (3,1,2),])

Header1 <- paste ("\\toprule \\toprule & \\multicolumn{2}{c}{\\textbf{\\emph{\\# lists seats}}} & \\multicolumn{2}{c}{\\textbf{\\emph{ENPS}}} & \\multicolumn{2}{c}{\\textbf{\\emph{Gallagher index}}} \\\\ \n")
Header2 <- paste ("[0.5ex] \\multicolumn{1}{l}{} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small}& \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} \\\\ \n")
Header3 <- paste ("\\multicolumn{1}{l}{\\textbf{(a) Main effect}} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(b) Heterogeneous effect}} \\\\ \\midrule \n")
Header5 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(c) Main effect (rank)}} \\\\ \\midrule \n")
Header6 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(b) Heterogeneous effect (rank)}} \\\\ \\midrule \n")

add.nobs <- paste ("[2.5ex] \\midrule num. obs & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,6],0)), "} \\\\ \n")
add.provs <- paste ("provinces & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,6],0)), "} \\\\ \n")
add.time <- paste ("elections & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,6],0)), "} \\\\ \n")

addtorow <- list ()
addtorow$pos <- list ()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 3
addtorow$pos[[5]] <- 7
addtorow$pos[[6]] <- 9
addtorow$pos[[7]] <- 13
addtorow$pos[[8]] <- 13
addtorow$pos[[9]] <- 13
addtorow$command <- c (Header1, Header2, Header3
                       , Header4, Header5, Header6
                       , add.nobs, add.provs, add.time)
print (xtable ( cbind (rows, tab01all)
                , align=c ("l","l","c","c","c","c","c","c")
                , digits=digits
                , caption="\\emph{District magnitude and the distribution of seats in Argentina, 1985-2015}"
                , label="T:results01")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="small"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )



### Replicating Table A4: Decomposition of the mechanical and psychological effect(s)
(tab03coef <- cbind (
  c (coef.m311a.fe[1,1], coef.m411a.fe[1,1], coef.m511a.fe[1,1], coef.m611a.fe[1,1], coef.m611a.fe[1,1]-coef.m311a.fe[1,1]
     , coef.m312a.fe[1:2,1], coef.m412a.fe[1:2,1], coef.m512a.fe[1:2,1], coef.m612a.fe[1:2,1], ((coef.m612a.fe[1,1]+coef.m612a.fe[2,1]*v.high)-(coef.m612a.fe[1,1]+coef.m612a.fe[2,1]*v.low)) - ((coef.m312a.fe[1,1]+coef.m312a.fe[2,1]*v.high)-(coef.m312a.fe[1,1]+coef.m312a.fe[2,1]*v.low)))
  , c (coef.m311d.fe[1,1], coef.m411d.fe[1,1], coef.m511d.fe[1,1], coef.m611d.fe[1,1], coef.m611d.fe[1,1]-coef.m311d.fe[1,1]
       , coef.m312d.fe[1:2,1], coef.m412d.fe[1:2,1], coef.m512d.fe[1:2,1], coef.m612d.fe[1:2,1], ((coef.m612d.fe[1,1]+coef.m612d.fe[2,1]*v.high)-(coef.m612d.fe[1,1]+coef.m612d.fe[2,1]*v.low)) - ((coef.m312d.fe[1,1]+coef.m312d.fe[2,1]*v.high)-(coef.m312d.fe[1,1]+coef.m312d.fe[2,1]*v.low)))
  
  , c (coef.m321a.fe[1,1], coef.m421a.fe[1,1], coef.m521a.fe[1,1], coef.m621a.fe[1,1], coef.m621a.fe[1,1]-coef.m321a.fe[1,1]
       , coef.m322a.fe[1:2,1], coef.m422a.fe[1:2,1], coef.m522a.fe[1:2,1], coef.m622a.fe[1:2,1], ((coef.m622a.fe[1,1]+coef.m622a.fe[2,1]*v.high)-(coef.m622a.fe[1,1]+coef.m622a.fe[2,1]*v.low)) - ((coef.m322a.fe[1,1]+coef.m322a.fe[2,1]*v.high)-(coef.m322a.fe[1,1]+coef.m322a.fe[2,1]*v.low)))
  , c (coef.m321d.fe[1,1], coef.m421d.fe[1,1], coef.m521d.fe[1,1], coef.m621d.fe[1,1], coef.m621d.fe[1,1]-coef.m321d.fe[1,1]
       , coef.m322d.fe[1:2,1], coef.m422d.fe[1:2,1], coef.m522d.fe[1:2,1], coef.m622d.fe[1:2,1], ((coef.m622d.fe[1,1]+coef.m622d.fe[2,1]*v.high)-(coef.m622d.fe[1,1]+coef.m622d.fe[2,1]*v.low)) - ((coef.m322d.fe[1,1]+coef.m322d.fe[2,1]*v.high)-(coef.m322d.fe[1,1]+coef.m322d.fe[2,1]*v.low)))
  
  , c (coef.m331a.fe[1,1], coef.m431a.fe[1,1], coef.m531a.fe[1,1], coef.m631a.fe[1,1], coef.m631a.fe[1,1]-coef.m331a.fe[1,1]
       , coef.m332a.fe[1:2,1], coef.m432a.fe[1:2,1], coef.m532a.fe[1:2,1], coef.m632a.fe[1:2,1], ((coef.m632a.fe[1,1]+coef.m632a.fe[2,1]*v.high)-(coef.m632a.fe[1,1]+coef.m632a.fe[2,1]*v.low)) - ((coef.m332a.fe[1,1]+coef.m332a.fe[2,1]*v.high)-(coef.m332a.fe[1,1]+coef.m332a.fe[2,1]*v.low)))
  , c (coef.m331d.fe[1,1], coef.m431d.fe[1,1], coef.m531d.fe[1,1], coef.m631d.fe[1,1], coef.m631d.fe[1,1]-coef.m331d.fe[1,1]
       , coef.m332d.fe[1:2,1], coef.m432d.fe[1:2,1], coef.m532d.fe[1:2,1], coef.m632d.fe[1:2,1], ((coef.m632d.fe[1,1]+coef.m632d.fe[2,1]*v.high)-(coef.m632d.fe[1,1]+coef.m632d.fe[2,1]*v.low)) - ((coef.m332d.fe[1,1]+coef.m332d.fe[2,1]*v.high)-(coef.m332d.fe[1,1]+coef.m332d.fe[2,1]*v.low)))
))
(tab03se <- cbind (
  c (coef.m311a.fe[1,2], coef.m411a.fe[1,2], coef.m511a.fe[1,2], coef.m611a.fe[1,2], NA
     , coef.m312a.fe[1:2,2], coef.m412a.fe[1:2,2], coef.m512a.fe[1:2,2], coef.m612a.fe[1:2,2], NA)
  , c (coef.m311d.fe[1,2], coef.m411d.fe[1,2], coef.m511d.fe[1,2], coef.m611d.fe[1,2], NA
       , coef.m312d.fe[1:2,2], coef.m412d.fe[1:2,2], coef.m512d.fe[1:2,2], coef.m612d.fe[1:2,2], NA)
  
  , c (coef.m321a.fe[1,2], coef.m421a.fe[1,2], coef.m521a.fe[1,2], coef.m621a.fe[1,2], NA
       , coef.m322a.fe[1:2,2], coef.m422a.fe[1:2,2], coef.m522a.fe[1:2,2], coef.m622a.fe[1:2,2], NA)
  , c (coef.m321d.fe[1,2], coef.m421d.fe[1,2], coef.m521d.fe[1,2], coef.m621d.fe[1,2], NA
       , coef.m322d.fe[1:2,2], coef.m422d.fe[1:2,2], coef.m522d.fe[1:2,2], coef.m622d.fe[1:2,2], NA)
  
  , c (coef.m331a.fe[1,2], coef.m431a.fe[1,2], coef.m531a.fe[1,2], coef.m631a.fe[1,2], NA
       , coef.m332a.fe[1:2,2], coef.m432a.fe[1:2,2], coef.m532a.fe[1:2,2], coef.m632a.fe[1:2,2], NA)
  , c (coef.m331d.fe[1,2], coef.m431d.fe[1,2], coef.m531d.fe[1,2], coef.m631d.fe[1,2], NA
       , coef.m332d.fe[1:2,2], coef.m432d.fe[1:2,2], coef.m532d.fe[1:2,2], coef.m632d.fe[1:2,2], NA)
))

# estimates with 95% CIs below
(tab03ci.full <- paste ("["
                        , sprintf ("%.2f", round (c (tab03coef[,c (1,3,5)]) + c (tab03se[,c (1,3,5)])*(-1)*tval, digits))
                        , ":"
                        , sprintf ("%.2f", round (c (tab03coef[,c (1,3,5)]) + c (tab03se[,c (1,3,5)])*(1)*tval, digits))
                        , "]", sep=""))
(tab03ci.sm <- paste ("["
                      , sprintf ("%.2f", round (c (tab03coef[,c (2,4,6)]) + c (tab03se[,c (2,4,6)])*(-1)*tval, digits))
                      , ":"
                      , sprintf ("%.2f", round (c (tab03coef[,c (2,4,6)]) + c (tab03se[,c (2,4,6)])*(1)*tval, digits))
                      , "]", sep=""))

# bootstrapped differences between coefficients
(tab.difs <- rbind (
  
  ## full sample
  
  # difference and 95% CIs $\\rightarrow$ main models
  cbind (
    c (coef.m611a.fe[1,1] - coef.m311a.fe[1,1], coef.m621a.fe[1,1] - coef.m321a.fe[1,1], coef.m631a.fe[1,1] - coef.m331a.fe[1,1])
    #, colMeans (coef.full[,c (1,7,13)] - coef.full[,c (2,8,14)])
    , apply (coef.full[,c (1,7,13)] - coef.full[,c (2,8,14)], 2, median)
    , apply (coef.full[,c (1,7,13)] - coef.full[,c (2,8,14)], 2, FUN=quantile, probs=.025)
    , apply (coef.full[,c (1,7,13)] - coef.full[,c (2,8,14)], 2, FUN=quantile, probs=.975)
  )
  
  # difference and 95% CIs $\\rightarrow$ heterogeneous effects
  , cbind (
    c (
      ((coef.m612a.fe[1,1]+coef.m612a.fe[2,1]*v.high) - (coef.m612a.fe[1,1]+coef.m612a.fe[2,1]*v.low)) - ((coef.m312a.fe[1,1]+coef.m312a.fe[2,1]*v.high) - (coef.m312a.fe[1,1]+coef.m312a.fe[2,1]*v.low))
      , ((coef.m622a.fe[1,1]+coef.m622a.fe[2,1]*v.high) - (coef.m622a.fe[1,1]+coef.m622a.fe[2,1]*v.low)) - ((coef.m322a.fe[1,1]+coef.m322a.fe[2,1]*v.high) - (coef.m322a.fe[1,1]+coef.m322a.fe[2,1]*v.low))
      , ((coef.m632a.fe[1,1]+coef.m632a.fe[2,1]*v.high) - (coef.m632a.fe[1,1]+coef.m632a.fe[2,1]*v.low)) - ((coef.m332a.fe[1,1]+coef.m332a.fe[2,1]*v.high) - (coef.m332a.fe[1,1]+coef.m332a.fe[2,1]*v.low)))
    , apply (
      ((coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.high) - (coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.low))
      -
        ((coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.high) - (coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.low))
      , 2, median)
    , apply (
      ((coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.high) - (coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.low))
      -
        ((coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.high) - (coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.low))
      , 2, FUN=quantile, probs=.025)
    , apply (
      ((coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.high) - (coef.full[,c (3,9,15)]+coef.full[,c (4,10,16)]*v.low))
      -
        ((coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.high) - (coef.full[,c (5,11,17)]+coef.full[,c (6,12,18)]*v.low))
      , 2, FUN=quantile, probs=.975))
  
  
  ## small-province sample
  
  # difference and 95% CIs $\\rightarrow$ main models
  , cbind (
    c (coef.m611d.fe[1,1] - coef.m311d.fe[1,1], coef.m621d.fe[1,1] - coef.m321d.fe[1,1], coef.m631d.fe[1,1] - coef.m331d.fe[1,1])
    , apply (coef.sm[,c (1,7,13)] - coef.sm[,c (2,8,14)], 2, median)
    , apply (coef.sm[,c (1,7,13)] - coef.sm[,c (2,8,14)], 2, FUN=quantile, probs=.025)
    , apply (coef.sm[,c (1,7,13)] - coef.sm[,c (2,8,14)], 2, FUN=quantile, probs=.975)
  )
  
  # difference and 95% CIs $\\rightarrow$ heterogeneous effects
  , cbind (
    c (
      ((coef.m612d.fe[1,1]+coef.m612d.fe[2,1]*v.high) - (coef.m612d.fe[1,1]+coef.m612d.fe[2,1]*v.low)) - ((coef.m312d.fe[1,1]+coef.m312d.fe[2,1]*v.high) - (coef.m312d.fe[1,1]+coef.m312d.fe[2,1]*v.low))
      , ((coef.m622d.fe[1,1]+coef.m622d.fe[2,1]*v.high) - (coef.m622d.fe[1,1]+coef.m622d.fe[2,1]*v.low)) - ((coef.m322d.fe[1,1]+coef.m322d.fe[2,1]*v.high) - (coef.m322d.fe[1,1]+coef.m322d.fe[2,1]*v.low))
      , ((coef.m632d.fe[1,1]+coef.m632d.fe[2,1]*v.high) - (coef.m632d.fe[1,1]+coef.m632d.fe[2,1]*v.low)) - ((coef.m332d.fe[1,1]+coef.m332d.fe[2,1]*v.high) - (coef.m332d.fe[1,1]+coef.m332d.fe[2,1]*v.low)))
    , apply (
      ((coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.high) - (coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.low))
      -
        ((coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.high) - (coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.low))
      , 2, median)
    , apply (
      ((coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.high) - (coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.low))
      -
        ((coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.high) - (coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.low))
      , 2, FUN=quantile, probs=.025)
    , apply (
      ((coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.high) - (coef.sm[,c (3,9,15)]+coef.sm[,c (4,10,16)]*v.low))
      -
        ((coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.high) - (coef.sm[,c (5,11,17)]+coef.sm[,c (6,12,18)]*v.low))
      , 2, FUN=quantile, probs=.975))
))
(tab.difs2 <- as.data.frame (cbind (
  c (rep ("full", 6), rep ("small", 6))
  , rep (c (rep ("main", 3), rep ("heterogeneous", 3)), 2)
)))
(tab.difs <- cbind (tab.difs2, tab.difs))
colnames (tab.difs) <- c ("sample", "effect", "dif", "dif.boot", "ci.low", "ci.high")

difs.ci <- with (tab.difs, paste ("[", sprintf ("%.2f", round (ci.low, 3)), ":", sprintf ("%.2f", round (ci.high, 3)), "]", sep=""))
tab03ci.full[tab03ci.full=="[NA:NA]"] <- difs.ci[c (1,4,2,5,3,6)]
tab03ci.sm[tab03ci.sm=="[NA:NA]"] <- difs.ci[c (1,4,2,5,3,6)+6]

# creating the table
odd <- seq (1, 28, by=2); even <- seq (2, 28, by=2)
(tab03all <- rbind (
  matrix (sprintf ("%.2f", round (tab03coef, digits)), nrow=14, byrow=F)
  , matrix (c (tab03ci.full, tab03ci.sm)[c (1:14, 43:56, 15:28, 57:70, 29:42, 71:84)], nrow=14, byrow=F)
)[sort (rep (1:14, 2)) + rep (c (0, 14), 14),])

rows <- c (
  # unconditional specifications
  "Mechanical effect: A $\\rightarrow$ B", ""
  , "Psychological effect (I): B $\\rightarrow$ D", ""
  , "Psychological effect (II): A $\\rightarrow$ C", ""
  , "Psychological effect (III): C $\\rightarrow$ D", ""
  , "Difference: [C $\\rightarrow$ D] - [A $\\rightarrow$ B]", ""
  
  # heterogeneous effects
  , "Mechanical effect: A $\\rightarrow$ B", "", "Mechanical effect: A $\\rightarrow$ B", "~~~~~$\\times$ \\emph{vote third party}"
  , "Psychological effect (I): B $\\rightarrow$ D", "", "Psychological effect (I): B $\\rightarrow$ D", "~~~~~$\\times$ \\emph{vote third party}"
  , "Psychological effect (II): A $\\rightarrow$ C", "", "Psychological effect (II): A $\\rightarrow$ C", "~~~~~$\\times$ \\emph{vote third party}"
  , "Psychological effect (III): C $\\rightarrow$ D", "", "Psychological effect (III): C $\\rightarrow$ D", "~~~~~$\\times$ \\emph{vote third party}"
  , "Difference: [C $\\rightarrow$ D] - [A $\\rightarrow$ B]", "")

(n.obs <- cbind (unlist (pdim (m311a.fe)$nT), unlist (pdim (m311d.fe)$nT), unlist (pdim (m321a.fe)$nT), unlist (pdim (m321d.fe)$nT), unlist (pdim (m331a.fe)$nT), unlist (pdim (m331d.fe)$nT))[c (3,1,2),])

Header1 <- paste ("\\toprule \\toprule & \\multicolumn{2}{c}{\\textbf{\\emph{\\# lists seats}}} & \\multicolumn{2}{c}{\\textbf{\\emph{ENPS}}} & \\multicolumn{2}{c}{\\textbf{\\emph{Gallagher index}}} \\\\ \n")
Header2 <- paste ("[0.5ex] \\multicolumn{1}{l}{} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} & \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small}& \\multicolumn{1}{c}{full} & \\multicolumn{1}{c}{small} \\\\ \n")
Header3 <- paste ("\\multicolumn{1}{l}{\\textbf{(a) Main effect}} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} & \\multicolumn{1}{c}{sample} & \\multicolumn{1}{c}{provinces} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{7}{l}{\\textbf{(b) Heterogeneous effect}} \\\\ \\midrule \n")

add.nobs <- paste ("[2.5ex] \\midrule num. obs & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[1,6],0)), "} \\\\ \n")
add.provs <- paste ("provinces & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[2,6],0)), "} \\\\ \n")
add.time <- paste ("elections & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs[3,6],0)), "} \\\\ \n")

addtorow <- list ()
addtorow$pos <- list ()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 10
addtorow$pos[[5]] <- 28
addtorow$pos[[6]] <- 28
addtorow$pos[[7]] <- 28
addtorow$command <- c (Header1, Header2, Header3
                       , Header4
                       , add.nobs, add.provs, add.time)
print (xtable ( cbind (rows, tab03all)
                , align=c ("l","l","c","c","c","c","c","c")
                , digits=digits
                , caption="\\emph{Contribution of the mechanical and psychological effects to the distribution of seats in Argentina, 1985-2015}"
                , label="T:results03")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="footnotesize"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )



### Replicating Table A5: Placebo tests
(tab04coef <- rbind (
  c (coef.m811a.fe[1,1], coef.m821a.fe[1,1], coef.m831a.fe[1,1], coef.m841a.fe[1,1], coef.m851a.fe[1,1], coef.m861a.fe[1,1], coef.m871a.fe[1,1], coef.m881a.fe[1,1])
  , c (coef.m811b.fe[1,1], coef.m821b.fe[1,1], coef.m831b.fe[1,1], coef.m841b.fe[1,1], coef.m851b.fe[1,1], coef.m861b.fe[1,1], coef.m871b.fe[1,1], coef.m881b.fe[1,1])
))
(tab04se <- rbind (
  c (coef.m811a.fe[1,2], coef.m821a.fe[1,2], coef.m831a.fe[1,2], coef.m841a.fe[1,2], coef.m851a.fe[1,2], coef.m861a.fe[1,2], coef.m871a.fe[1,2], coef.m881a.fe[1,2])
  , c (coef.m811b.fe[1,2], coef.m821b.fe[1,2], coef.m831b.fe[1,2], coef.m841b.fe[1,2], coef.m851b.fe[1,2], coef.m861b.fe[1,2], coef.m871b.fe[1,2], coef.m881b.fe[1,2])
))

# estimates with 95% CIs below
(tab04ci.full <- paste ("["
                        , sprintf ("%.2f", round (c (tab04coef) + c (tab04se)*(-1)*tval, digits))
                        , ":"
                        , sprintf ("%.2f", round (c (tab04coef) + c (tab04se)*(1)*tval, digits))
                        , "]", sep=""))

# creating the table
odd <- seq (1, 4, by=2); even <- seq (2, 4, by=2)
(tab04all <- rbind (
  matrix (sprintf ("%.2f", round (tab04coef, digits)), nrow=2, byrow=F)
  , matrix (c (tab04ci.full), nrow=2, byrow=F))[c (1,3,2,4),])

rows <- c ("\\emph{magnitude}", "", "\\emph{magnitude}", "")

(n.obs01 <- cbind (unlist (pdim (m811a.fe)$nT), unlist (pdim (m821a.fe)$nT), unlist (pdim (m831a.fe)$nT), unlist (pdim (m841a.fe)$nT), unlist (pdim (m851a.fe)$nT), unlist (pdim (m861a.fe)$nT), unlist (pdim (m871a.fe)$nT), unlist (pdim (m881a.fe)$nT))[c (3,1,2),])
(n.obs02 <- cbind (unlist (pdim (m811b.fe)$nT), unlist (pdim (m821b.fe)$nT), unlist (pdim (m831b.fe)$nT), unlist (pdim (m841b.fe)$nT), unlist (pdim (m851b.fe)$nT), unlist (pdim (m861b.fe)$nT), unlist (pdim (m871b.fe)$nT), unlist (pdim (m881b.fe)$nT))[c (3,1,2),])

Header1 <- paste ("\\toprule \\toprule & \\multicolumn{1}{c}{\\textbf{\\emph{revenues}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{public}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{infant}}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{1}{c}{\\textbf{\\emph{per capita}}} & \\multicolumn{1}{c}{\\textbf{\\emph{\\% own}}} & \\multicolumn{1}{c}{\\textbf{\\emph{}}} & \\multicolumn{1}{c}{\\textbf{\\emph{\\% automatic}}} & \\multicolumn{1}{c}{\\textbf{\\emph{\\% discretionary}}} & \\multicolumn{1}{c}{\\textbf{\\emph{employees}}} & \\multicolumn{1}{c}{\\textbf{\\emph{unemployment}}} & \\multicolumn{1}{c}{\\textbf{\\emph{mortality}}} \\\\ \n")
Header3 <- paste ("\\multicolumn{1}{l}{\\textbf{(a) Full sample}} & \\multicolumn{1}{c}{\\textbf{\\emph{(log)}}} & \\multicolumn{1}{c}{\\textbf{\\emph{revenues}}} & \\multicolumn{1}{c}{\\textbf{\\emph{\\% royalties}}} & \\multicolumn{1}{c}{\\textbf{\\emph{transfers}}} & \\multicolumn{1}{c}{\\textbf{\\emph{transfers}}} & \\multicolumn{1}{c}{\\textbf{\\emph{(per 1,000)}}} & \\multicolumn{1}{c}{\\textbf{\\emph{rate (\\%)}}} & \\multicolumn{1}{c}{\\textbf{\\emph{(per 1,000)}}} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{9}{l}{\\textbf{(b) Small provinces}} \\\\ \\midrule \n")

add.nobs01 <- paste ("[2.5ex] \\midrule num. obs & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[1,8],0)), "} \\\\ \n")
add.provs01 <- paste ("provinces & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[2,8],0)), "} \\\\ \n")
add.time01 <- paste ("elections & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs01[3,8],0)), "} \\\\ \n")
add.nobs02 <- paste ("[2.5ex] \\midrule num. obs & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[1,8],0)), "} \\\\ \n")
add.provs02 <- paste ("provinces & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[2,8],0)), "} \\\\ \n")
add.time02 <- paste ("elections & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,1],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,2],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,3],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,4],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,5],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,6],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,7],0)), "} & \\multicolumn{1}{c}{", sprintf ("%.0f", round (n.obs02[3,8],0)), "} \\\\ \n")

addtorow <- list ()
addtorow$pos <- list ()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 2
addtorow$pos[[5]] <- 2
addtorow$pos[[6]] <- 2
addtorow$pos[[7]] <- 2
addtorow$pos[[8]] <- 4
addtorow$pos[[9]] <- 4
addtorow$pos[[10]] <- 4
addtorow$command <- c (Header1, Header2, Header3
                       , add.nobs01, add.provs01, add.time01
                       , Header4
                       , add.nobs02, add.provs02, add.time02)
print (xtable ( cbind (rows, tab04all)
                , align=c ("l","l","c","c","c","c","c","c","c","c")
                , digits=digits
                , caption="\\emph{Placebo tests. The effect of district magnitude on time-varying pseudo-outcomes in Argentina, 1985-2011}"
                , label="T:results04")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="small"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )





########## BUILDING THE FIGURES

## Setting values that we will keep constant when presenting the results
cex.points <- 3.5
cex.text <- 21
col.line <- "gray66"
col.soft <- "gray60"
cex.bars <- 0.045*3
pd <- position_dodge (width = 0.5)
lim.x <- c (0,1)

### Calculating the values reported in Figure 1 (mchanical vs. psychological effect)

## Larger magnitudes in concurrent years
sort (unique (baseFull[!is.na (baseFull$largeMidterm) & baseFull$largeMidterm==0,]$district)) ## CHU  CTES FSA  LR   MIS  RN   SFE  SGO  STA  TDF  TUC

## group A
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="A",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="A",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="A",]$dispGall), 2)

## group B
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="B",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="B",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="B",]$dispGall), 2)

## group C
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="C",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="C",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="C",]$dispGall), 2)

## group D
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="D",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="D",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==0 & baseFull$ffCode=="D",]$dispGall), 2)

## Larger magnitudes in midterm years
sort (unique (baseFull[!is.na (baseFull$largeMidterm) & baseFull$largeMidterm==1,]$district)) ## CABA CAT  CHA  ER   LP   NQN  SC   SL

## group A
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="A",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="A",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="A",]$dispGall), 2)

## group B
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="B",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="B",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="B",]$dispGall), 2)

## group C
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="C",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="C",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="C",]$dispGall), 2)

## group D
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="D",]$nListSeats), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="D",]$enps), 2)
round (mean (baseFull[!is.na (baseFull$ffCode) & baseFull$odd==1 & baseFull$largeMidterm==1 & baseFull$ffCode=="D",]$dispGall), 2)


### Replicating Figure 3: Electoral coordination
(dplot02a <- as.data.frame (cbind (
  
  # estimates
  c (tab02coef[1,], (tab02coef[2,]+tab02coef[3,]*v.low), (tab02coef[2,]+tab02coef[3,]*v.high))
  
  # SEs
  , c (tab02se[1,], sqrt (1^2*vcov.m212a.fe[1,1] + v.low^2*vcov.m212a.fe[2,2] + 2*v.low*vcov.m212a.fe[1,2])
       , sqrt (1^2*vcov.m212d.fe[1,1] + v.low^2*vcov.m212d.fe[2,2] + 2*v.low*vcov.m212d.fe[1,2])
       , sqrt (1^2*vcov.m222a.fe[1,1] + v.low^2*vcov.m222a.fe[2,2] + 2*v.low*vcov.m222a.fe[1,2])
       , sqrt (1^2*vcov.m222d.fe[1,1] + v.low^2*vcov.m222d.fe[2,2] + 2*v.low*vcov.m222d.fe[1,2])
       , sqrt (1^2*vcov.m232a.fe[1,1] + v.low^2*vcov.m232a.fe[2,2] + 2*v.low*vcov.m232a.fe[1,2])
       , sqrt (1^2*vcov.m232d.fe[1,1] + v.low^2*vcov.m232d.fe[2,2] + 2*v.low*vcov.m232d.fe[1,2])
       , sqrt (1^2*vcov.m212a.fe[1,1] + v.high^2*vcov.m212a.fe[2,2] + 2*v.high*vcov.m212a.fe[1,2])
       , sqrt (1^2*vcov.m212d.fe[1,1] + v.high^2*vcov.m212d.fe[2,2] + 2*v.high*vcov.m212d.fe[1,2])
       , sqrt (1^2*vcov.m222a.fe[1,1] + v.high^2*vcov.m222a.fe[2,2] + 2*v.high*vcov.m222a.fe[1,2])
       , sqrt (1^2*vcov.m222d.fe[1,1] + v.high^2*vcov.m222d.fe[2,2] + 2*v.high*vcov.m222d.fe[1,2])
       , sqrt (1^2*vcov.m232a.fe[1,1] + v.high^2*vcov.m232a.fe[2,2] + 2*v.high*vcov.m232a.fe[1,2])
       , sqrt (1^2*vcov.m232d.fe[1,1] + v.high^2*vcov.m232d.fe[2,2] + 2*v.high*vcov.m232d.fe[1,2])
  ))))
colnames (dplot02a) <- c ("estimate", "se")
dplot02a$tval <- NA
dplot02a[seq (1, nrow (dplot02a), by=2),]$tval <- tval
dplot02a[seq (2, nrow (dplot02a), by=2),]$tval <- tval.sm
dplot02a$low95 <- with (dplot02a, estimate - tval*se)
dplot02a$high95 <- with (dplot02a, estimate + tval*se)
dplot02a$sample <- factor (rep (c ("full sample", "small provinces"), nrow (dplot02a)/2))
dplot02a$dv <- factor (rep (c (rep ("(a) # lists running", 2), rep ("(b) ENPV", 2), rep ("(c) vote first two", 2)), 3))
dplot02a$dv <- factor (dplot02a$dv, levels=c ("(a) # lists running", "(b) ENPV", "(c) vote first two"))
dplot02a$effect <- factor (c (rep ("main\neffect", 6), rep ("vote third\nparty: 5%", 6), rep ("vote third\nparty: 15%", 6)))
dplot02a$effect <- factor (dplot02a$effect, levels=c ("main\neffect", "vote third\nparty: 5%", "vote third\nparty: 15%"))

(plot02a <- ggplot (dplot02a, aes (y=estimate, x=effect, colour=sample, group=sample, fill=sample, shape=sample))
 + geom_hline (yintercept=0, colour=col.line, linetype=2)
 + geom_point (size=cex.points, position=pd)
 + geom_errorbar (aes (ymin=low95, ymax=high95), position=pd, size=cex.points/4, width=cex.bars)
 + scale_colour_grey () ## grayscale for the printed version; comment out for colored version
 + ylab ("Marginal effect of magnitude")
 + xlab ("")
 + theme_bw ()
 + facet_wrap ( ~ dv, scales="free_y")
 + theme (text=element_text (size=cex.text), strip.text.y=element_text (size=cex.text), strip.text.x=element_text (size=cex.text))
 + theme (legend.position = "bottom", legend.title = element_blank ())
 )



### Replicating Figure 4: Aggregate outcomes
(dplot01a <- as.data.frame (cbind (
  
  # estimates
  c (tab01coef[1,], (tab01coef[2,]+tab01coef[3,]*v.low), (tab01coef[2,]+tab01coef[3,]*v.high))
  
  # SEs
  , c (tab01se[1,], sqrt (1^2*vcov.m112a.fe[1,1] + v.low^2*vcov.m112a.fe[2,2] + 2*v.low*vcov.m112a.fe[1,2])
       , sqrt (1^2*vcov.m112d.fe[1,1] + v.low^2*vcov.m112d.fe[2,2] + 2*v.low*vcov.m112d.fe[1,2])
       , sqrt (1^2*vcov.m122a.fe[1,1] + v.low^2*vcov.m122a.fe[2,2] + 2*v.low*vcov.m122a.fe[1,2])
       , sqrt (1^2*vcov.m122d.fe[1,1] + v.low^2*vcov.m122d.fe[2,2] + 2*v.low*vcov.m122d.fe[1,2])
       , sqrt (1^2*vcov.m132a.fe[1,1] + v.low^2*vcov.m132a.fe[2,2] + 2*v.low*vcov.m132a.fe[1,2])
       , sqrt (1^2*vcov.m132d.fe[1,1] + v.low^2*vcov.m132d.fe[2,2] + 2*v.low*vcov.m132d.fe[1,2])
       , sqrt (1^2*vcov.m112a.fe[1,1] + v.high^2*vcov.m112a.fe[2,2] + 2*v.high*vcov.m112a.fe[1,2])
       , sqrt (1^2*vcov.m112d.fe[1,1] + v.high^2*vcov.m112d.fe[2,2] + 2*v.high*vcov.m112d.fe[1,2])
       , sqrt (1^2*vcov.m122a.fe[1,1] + v.high^2*vcov.m122a.fe[2,2] + 2*v.high*vcov.m122a.fe[1,2])
       , sqrt (1^2*vcov.m122d.fe[1,1] + v.high^2*vcov.m122d.fe[2,2] + 2*v.high*vcov.m122d.fe[1,2])
       , sqrt (1^2*vcov.m132a.fe[1,1] + v.high^2*vcov.m132a.fe[2,2] + 2*v.high*vcov.m132a.fe[1,2])
       , sqrt (1^2*vcov.m132d.fe[1,1] + v.high^2*vcov.m132d.fe[2,2] + 2*v.high*vcov.m132d.fe[1,2])
  ))))
colnames (dplot01a) <- c ("estimate", "se")
dplot01a$tval <- NA
dplot01a[seq (1, nrow (dplot01a), by=2),]$tval <- tval
dplot01a[seq (2, nrow (dplot01a), by=2),]$tval <- tval.sm
dplot01a$low95 <- with (dplot01a, estimate - tval*se)
dplot01a$high95 <- with (dplot01a, estimate + tval*se)
dplot01a$sample <- factor (rep (c ("full sample", "small provinces"), nrow (dplot01a)/2))
dplot01a$dv <- factor (rep (c (rep ("(a) # lists seats", 2), rep ("(b) ENPS", 2), rep ("(c) Gallagher index", 2)), 3))
dplot01a$dv <- factor (dplot01a$dv, levels=c ("(a) # lists seats", "(b) ENPS", "(c) Gallagher index"))
dplot01a$effect <- factor (c (rep ("main\neffect", 6), rep ("vote third\nparty: 5%", 6), rep ("vote third\nparty: 15%", 6)))
dplot01a$effect <- factor (dplot01a$effect, levels=c ("main\neffect", "vote third\nparty: 5%", "vote third\nparty: 15%"))

(plot01a <- ggplot (dplot01a, aes (y=estimate, x=effect, colour=sample, group=sample, fill=sample, shape=sample))
 + geom_hline (yintercept=0, colour=col.line, linetype=2)
 + geom_point (size=cex.points, position=pd)
 + geom_errorbar (aes (ymin=low95, ymax=high95), position=pd, size=cex.points/4, width=cex.bars)
 + scale_colour_grey () ## grayscale for the printed version; comment out for colored version
 + ylab ("Marginal effect of magnitude")
 + xlab ("")
 + theme_bw ()
 + facet_wrap ( ~ dv, scales="free_y")
 + theme (text=element_text (size=cex.text), strip.text.y=element_text (size=cex.text), strip.text.x=element_text (size=cex.text))
 + theme (legend.position = "bottom", legend.title = element_blank ())
)



### Replicating Figure 5: Decomposing the mechanical and psychological effects
(dplot03a <- as.data.frame (cbind (
  
  # estimates
  c (tab03coef[1,], tab03coef[2,], tab03coef[3,], tab03coef[5,])
  
  # SEs
  , c (tab03se[1,], tab03se[2,], tab03se[3,], rep (NA, 6))
)))
colnames (dplot03a) <- c ("estimate", "se")
dplot03a$tval <- NA
dplot03a[seq (1, nrow (dplot03a), by=2),]$tval <- tval
dplot03a[seq (2, nrow (dplot03a), by=2),]$tval <- tval.sm
dplot03a$low95 <- with (dplot03a, estimate - tval*se)
dplot03a$high95 <- with (dplot03a, estimate + tval*se)
dplot03a[is.na (dplot03a$se),]$low95 <- tab.difs[c (1,7,2,8,3,9),]$ci.low
dplot03a[is.na (dplot03a$se),]$high95 <- tab.difs[c (1,7,2,8,3,9),]$ci.high
dplot03a$sample <- factor (rep (c ("full\nsample", "small\nprovinces"), nrow (dplot03a)/2))
# dplot03a$decomposition <- factor (c (rep ("Mechanical\neffect: A \u27F6 B", 6), rep ("Psychological\neffect (I): B \u27F6 D", 6), rep ("Psychological\neffect (II): A \u27F6 C", 6), rep ("Difference:\n[C \u27F6 D] - [A \u27F6 B]", 6)))
dplot03a$decomposition <- factor (c (rep ("Mechanical\neffect: A -> B", 6), rep ("Psychological\neffect (I): B -> D", 6), rep ("Psychological\neffect (II): A -> C", 6), rep ("Difference:\n[C -> D] - [A -> B]", 6)))
# dplot03a$decomposition <- factor (dplot03a$decomposition, levels=c ("Mechanical\neffect: A \u27F6 B", "Psychological\neffect (I): B \u27F6 D", "Psychological\neffect (II): A \u27F6 C", "Difference:\n[C \u27F6 D] - [A \u27F6 B]"))
dplot03a$decomposition <- factor (dplot03a$decomposition, levels=c ("Mechanical\neffect: A -> B", "Psychological\neffect (I): B -> D", "Psychological\neffect (II): A -> C", "Difference:\n[C -> D] - [A -> B]"))
dplot03a$dv <- factor (rep (c (rep ("(a) # lists seats", 2), rep ("(b) ENPS", 2), rep ("(c) Gallagher index", 2)), 4))
dplot03a$dv <- factor (dplot03a$dv, levels=c ("(a) # lists seats", "(b) ENPS", "(c) Gallagher index"))
dplot03a$effect <- factor ("main\neffect")

(plot03a <- ggplot (dplot03a, aes (y=estimate, x=sample, colour=sample, group=sample, fill=sample, shape=sample))
 + geom_hline (yintercept=0, colour=col.line, linetype=2)
 + geom_point (size=cex.points, position=pd)
 + geom_errorbar (aes (ymin=low95, ymax=high95), position=pd, size=cex.points/4, width=cex.bars)
 + scale_colour_grey () ## grayscale for the printed version; comment out for colored version
 + ylab ("Marginal effect of magnitude")
 + xlab ("")
 + theme_bw ()
 + facet_grid (dv ~ decomposition, scales="free_y")
 + theme (text=element_text (size=cex.text), strip.text.y=element_text (size=cex.text), strip.text.x=element_text (size=cex.text))
 + theme (legend.position = "none")
)



### Replicating Figure 6: Placebo tests
(dplot04a <- as.data.frame (cbind (
  c (tab04coef[1,], tab04coef[2,])  # estimates
  , c (tab04se[1,], tab04se[2,])  # SEs
)))
colnames (dplot04a) <- c ("estimate", "se")
dplot04a$tval <- NA
dplot04a[1:nrow (dplot04a)/2,]$tval <- tval
dplot04a[(nrow (dplot04a)/2+1):nrow (dplot04a),]$tval <- tval.sm
dplot04a$low95 <- with (dplot04a, estimate - tval*se)
dplot04a$high95 <- with (dplot04a, estimate + tval*se)
dplot04a$sample <- factor (c (rep ("full\nsample", nrow (dplot04a)/2), rep ("small\nprovinces", nrow (dplot04a)/2)))
dplot04a$dv <- factor (rep (c ("revenues\nper capita (log)", "% own\nrevenues", "% royalties", "% automatic\ntransfers", "% discretionary\ntransfers", "public employees\n(per 1,000)", "unemployment\nrate (%)", "infant mortality\n(per 1,000)"), 2))
dplot04a$dv <- factor (dplot04a$dv, levels=c ("revenues\nper capita (log)", "% own\nrevenues", "% royalties", "% automatic\ntransfers", "% discretionary\ntransfers", "public employees\n(per 1,000)", "unemployment\nrate (%)", "infant mortality\n(per 1,000)"))

(plot04a <- ggplot (dplot04a, aes (y=estimate, x=sample, colour=sample, group=sample, fill=sample, shape=sample))
 + geom_hline (yintercept=0, colour=col.line, linetype=2)
 + geom_point (size=cex.points, position=pd)
 + geom_errorbar (aes (ymin=low95, ymax=high95), position=pd, size=cex.points/4, width=cex.bars)
 + scale_colour_grey () ## grayscale for the printed version; comment out for colored version
 + ylab ("Marginal effect of magnitude")
 + xlab ("")
 + theme_bw ()
 + facet_wrap ( ~ dv, scales="free_y", nrow=2)
 + theme (text=element_text (size=cex.text), strip.text.y=element_text (size=cex.text), strip.text.x=element_text (size=cex.text))
 + theme (legend.position = "none")
 )


## Exporting the plots

# different heights for different numbers of panels:
plotWidth <- 1100
plotHeight3 <- (240*3)+50
plotHeight2 <- (240*2)+50
plotHeight1 <- (320*1)+50

png ("fig_results01.png", w=plotWidth, h=plotHeight1)
plot01a
dev.off ()

png ("fig_results02.png", w=plotWidth, h=plotHeight1)
plot02a
dev.off ()

png ("fig_results03.png", w=plotWidth, h=plotHeight3)
plot03a
dev.off ()

png ("fig_results04.png", w=plotWidth, h=plotHeight2)
plot04a
dev.off ()
